Accuracy of carotid strain estimates from ultrasonic wall tracking: a multiphysics model study

Introduction Ultrasound imaging of the carotid artery is a common procedure when screening for cardiovascular disease, as the vessel is particularly prone to atherosclerosis and easily accessible with ultrasound probes. However, to reveal abnormalities in carotid flow and wall deformation with increased sensitivity and specificity, improved imaging modalities are desired. Since in-vitro and in-vivo testing of new imaging algorithms often lack ground truth information, ultrasonic image simulation based on fully known and realistic vascular behavior can be expected to support image development. In this context, we developed a multiphysical simulation environment which integrates advanced numerical methods to calculate complex flow patterns and mechanical deformations on the one hand with an ultrasonic simulator on the other hand. For the ultrasonic image modeling, we relied on the Field II-software [ref], which can simulate with great scanning flexibility images of arbitrary tissue (i.e. both blood and arterial wall) by representing it as an ensemble of point scatterers on which the ultrasound waves reflect. In a first phase, we coupled computational fluid dynamics (CFD) with Field II, allowing simulation of radiofrequency (RF) data from realistic and complex flow fields, by moving the scatterers according to the spatially and temporally interpolated velocity fields obtained from CFD. We demonstrated the realism of the simulation environment with color flow imaging and pulsed wave Doppler examples in the carotid artery [ref] and by validating the simulation strategy on an in-vitro flow phantom of the carotid bifurcation [ref]. Fundamental limitations of this approach were the absence of the vessel wall signal and the rigid vessel walls. In a next phase, we integrated fluid-structure interaction (FSI) simulations with the ultrasound simulator [ref], which allowed to simultaneously assess the complex flow field and vessel wall deformation by coupling the numerical solution of a dedicated flow solver and structural solver. Hence, an FSI-Field II integration offers the possibility to simultaneously simulate the RF-signal of the blood pool and the moving vessel wall. This, however, requires a much more complex coupling methodology (derivation of the scatterer displacement) compared to the CFD-Field II coupling, due to the temporally varying fluid volume and the layered vessel wall, as explained in [ref]. We demonstrated the complex FSI-Field II coupling in a 3D straight tube, representative of the common carotid artery [ref]. In this paper, we further advance the realism of the synthetic vascular imaging set-up with the simulation of the 3D blood flow and arterial mechanics of a patient-specific carotid bifurcation model. The model includes the tissue surrounding the vascular wall, which is a stabilizing factor for the fluid-structure interaction simulation but also results in a more realistic ultrasonic visualization of the arterial territory. The extensive methodology behind the FSI-Field II coupling will be briefly discussed, before demonstrating the realism of our multiphysics modeling with simulated duplex images. Subsequently, we assess the accuracy of ultrasonic estimation of radial and circumferential strain and related material properties as derived from a wall tracking algorithm [ref], which is the main focus of this paper. The performance of wall tracking algorithms was previously analyzed in a simplified tube configuration (deforming in an axial symmetrical way) in order to ease the interpretation of the measured distensions (∆D) and the corresponding circumferential strain (∆D/D) estimates [ref]. We demonstrated that distension and circumferential strain estimates as obtained from wall tracking measurements are affected by the physics of the ultrasonic image formation and should be interpreted with care when linking them to the mechanical properties of the wall tissue [ref]. As wall tracking is applied to more complexly deforming blood vessels like the carotid artery, we anticipate that its intricate wall mechanics will even further complicate the analysis. Therefore, simulated distension and strain estimates (circumferential and radial components) derived from an ultrasonic wall tracking algorithm will be compared with the true mechanical deformation and complex 3D strains as known from the FSI-simulation. Methods 1. FSI-simulations 1.a Numerical approach A partitioned FSI-approach was followed, i.e. the equations for the flow and structural domain were solved separately with a dedicated flow (Fluent 12.0.16, Ansys, Canonsburg, PA, USA) and structural solver (Abaqus 6.7, Simulia,Inc., Providence, RI, USA). More information on the applied numerical approaches in these dedicated solvers can be found in [ref]. The solutions for the fluid and structural domain were coupled using in-house code (Tango) with Dirichlet-Neumann partitioning (the flow problem is solved for a given displacement of the fluid-structure interface, and the structural problem is solved for a given stress distribution on the wet side of the structure). A converged solution for both the fluid and structural equations and the coupling conditions was found by performing coupling iterations between both solvers, until equilibrium between the fluid and structure was reached. To enhance convergence of the coupling iterations, an Interface Quasi-Newton (IQN) method was used [ref]. Note that the fluid and structural domain have inherently different grid formulations, which was solved by using an Arbitrary Lagrangian Eulerian (ALE) method for the fluid domain. For more information on the IQN and ALE method, we refer to [ref]. 1.b (Meshing the) carotid geometry The 3D geometry was reconstructed from MRI-scans of a stenosed carotid bifurcation of an 83-year old volunteer. The MRI-scan sequence covered a 2 cm-region centered around the bifurcation. The geometry was artificially prolonged at the inand outlets (total length of the model was 6.5 cm), in order to obtain fully developed flow in the region of interest (the bifurcation), but also to create a sufficiently long computational phantom for the simulated ultrasonic scanning sequence. An in-house open-source mesh generation code (Pyformex) was used to construct a computational grid for the vascular wall. This allowed to create a 4-layered mesh of the vascular wall (to some extent mimicking the intima-media-adventitia layers), consisting of 31680 first order hexahedral elements. Pyformex was subsequently used to mesh the fluid domain with hexahedrons, resulting in 87522 elements and a matching grid at the fluid-structure interface. The structural domain was further expanded to also include the tissue surrounding the vascular wall. Although the model becomes more computationally expensive, adding the external tissue had a stabilizing effect for the FSI-simulations and allowed for more realistic ultrasound simulations. The carotid artery was embedded in a cylinder with a radius of 2 cm. Obviously, the surrounding tissue requires meshing as well, which is a challenging undertaking. Meshing this tissue with hexahedral elements is not an option since this element type causes intersections close to the bifurcation. As such, although hexahedrons are preferred for structural analysis (tetrahedrons behave stiffer), a first-order tetrahedral element type was chosen for the tissue volume. To ensure continuity of the hexahedral and tetrahedral meshes at the wall-tissue interface (to avoid interaction problems in the structural solver), the quadrilaterals at the outer surface of the vascular wall were split into triangles to match with the tetrahedrons of the tissue domain. A cross-section of the complete computational grid is shown in fig.1. 1.c FSI-setup Fluid domain: At the inand outlets of the carotid geometry, we imposed physiologically realistic boundary conditions. We measured a velocity profile with pulsed wave Doppler (12L linear array vascular probe, GE Medical Systems, Milwaukee, WI, USA) in the common carotid of a healthy volunteer, which was further applied as a mass flow inlet condition. Outflow percentages were imposed at the outlets (35% at the external and 65 % at the internal carotid). Since the absolute pressure level is undetermined for such a setup, we added a non-invasively measured pressure (varying in time, with a pulse pressure of 40 mmHg) to the obtained fluid pressure distribution. As such, when transferring the interface stress to the structural solver, a realistic pressure value was imposed on the wet side of the structure. Blood was modeled as a Newtonian liquid with a viscosity of 3.5 mPas and a density of 1050 kg/m3. Solid domain: Assuming that the mechanical properties of the vessel wall material can be linearized around the operating pressure, we modeled the vessel wall as a linear elastic material with Young’s modulus of 250 kPa, density of 1200 kg/m3 and Poisson modulus of 0.49 (nearly incompressible). The properties of the surrounding tissue were chosen to obtain a realistic distension degree for the chosen vessel elasticity: a Young’s modulus of 10 kPa and Poisson modulus of 0.3 [ref]. Longitudinal movement of the inand outlets was prevented. We refer to fig.1 for a complete overview of the applied fluid and solid boundary conditions. The cardiac cycle of 1s was divided into timesteps of 5 ms and 2 cycles were computed to obtain results independent of transient effects. The coupling algorithm was executed on one core, the flow solver on eight cores and the structural solver on eight cores of a dedicated machine with two Intel Xeon 2x Quad-core Intel Xeon processors (2.66 GHz). 2. Ultrasound-simulations 2.a Field II The Field II software [ref] was used to simulate the RF-signals from the fluid and structural domain. This modeling approach allows simulating arbitrary transducers and scansequencing with great flexibility and is based on the spatial impulse response estimation as described by Tupholme and Stepanishen [ref]. This simulation strategy is limited to linear wave propagation and determines the ultrasound field based on the ultrasonic excitation pulse, the temporal impulse responses of the transmitting and receiving transducers, and the spatial impulse response at a given point. For further details on the theoretical background, we refer to [ref]. The RFsignals can be simulated with a high degree of realism because Field II models tissue as a distribution of (random) point scatterers, whose position can be updated for each simulated ultrasound beam. By moving the scatterers according to flow fields and wall deformations obtained from CFD/FSI, imaging algorithms can be studied in complex conditions [refs]. Fluid phantom: We refer to [] for details on how the scatterers can be propagated using FSI simulation results accounting for the temporally varying fluid volume and the applied ALE grid formulation. [ref]. Wall phantom: Deriving the scatterer displacements for the wall phantom was less complex due to the Langrangian grid formulation, i.e. the grid displacement corresponds to the material (scatterer) displacement. However, to account for local changes in material (acoustic) properties, random point scatterers were generated for each element of the wall mesh. Further, as can be seen in echo images, the vessel wall also causes specular reflections due to transitions between different tissue types. These cannot be simulated but only mimicked in Field II, by placing scatterers in a structured fashion at the borders of the vessel wall (i.e. tissue/vessel wall and vessel wall/blood). For more details on the coupling methodology for the wall, we again refer to [ref]. Tissue phantom: As explained above, the model also includes the surrounding tissue. To reduce computational times, scatterers were not generated for each mesh element of the tissue (as for the wall phantom), but for the complete cylinder surrounding the carotid artery, comprising the fluid, wall and tissue domain. Afterwards, scatterers created inside the arterial wall and fluid volume were removed. Although the number of scatterers in Field II is determined by the resolution of the imaging system (10 scatterers per resolution cell assuring Gaussian distributed RF-signals), the total amount of tissue scatterers was reduced by a factor 10 for calculation purposes. This was justified by the fact that no visual differences were apparent between images with the full and reduced amount of scatterers. Further, no in-depth evaluation of tissue RF-signals was included in this study. 2.b Imaging setup Both for the duplex scanning and the wall tracking, a linear array transducer was modeled, with a focal depth position at 2 cm. Each transducer element was divided into four smaller rectangular mathematical elements so that the backscattered signal from each point scatterer was simulated with sufficient accuracy. A dynamic focus and expanding aperture was used on receive to retain constant imaging properties throughout depth. To reduce beam sidelobes, apodization was applied. Duplex scanning: A duplex scan is the superposition of a color flow image (CFI) on a B-mode image. However, the image acquisition requirements of CFI and B-mode are inherently different due to dissimilar spatial and temporal resolution requirements, and therefore compromises have to be made. A 5 MHz centre frequency was chosen for both image acquisitions, but the beam density was doubled for the B-mode imaging and the pulse length was increased for CFI from 1.5 to 4 pulse periods. To achieve these differing imaging properties, an interleaved scanning scheme was applied, switching between color flow and B-mode acquisitions. This resulted in a frame rate of 12 fps. Full details on the imaging setup can be found in table 1. The color flow velocity estimates were estimated with the autocorrelation method for phase-shift estimation, as proposed for ultrasound applications by Kasai et al [ref]. The axial velocity vz (cfr. coordinate system on fig.1) was calculated according to: ) )) 1 ( ( Re )) 1 ( Im( arctan( 4 ^ ^ 0 R R f cPRF vz   , with PRF (=pulse repetition frequency) the frequency of emitting ultrasound beams, f0 the centre frequency of the ultrasound pulse and ) 1 ( ^ R the estimated autocorrelation function at lag 1. The ) 1 ( ^ R estimate was averaged over an ensemble of 10 slow-time samples. The phantom was angled 70 degrees towards the ultrasound scanline to reduce flow transversal to the beam. To improve the CFI frame rate, we applied a beam interleaved acquisition scheme as described in [ref], typically used when the Doppler PRF (as determined by the imaged velocity range) is chosen lower than the maximal possible PRFmax (as determined by the imaging depth). We chose a setup with a PRF of 4 kHz, a PRFmax of 16 kHz, resulting in an interleave groupsize of 4 beams. Wall tracking: Vessel wall velocities were estimated with a modified autocorrelation approach [ref], allowing to determine vessel wall motion as: z[t+∆t]=z[t]+v[t]∆t, with z[t] the position in the vessel wall, v[t] the estimated velocity, and ∆t the velocity resolution corresponding to the packet size times the pulse repetition period (3•1/1000). A relatively low PRF of 1 kHz was chosen, because tissue velocities in the carotid artery are typically in a lower range than blood velocities. The phantom was not angled since, for this application, ultrasound beams are ideally emitted perpendicularly to the vessel wall. Compared to the duplex scanning, an imaging setup with a much higher resolution was required and the 12L linear array probe (GE Medical Systems, Milwaukee, WI, USA), as used in the applied distension software [ref] was modeled with a 1.5 period sinusoidal pulse excitation of 8 MHz centre frequency. The RF-signal from the surrounding tissue and fluid domain was neglected for this application. Further details on the imaging setup can be found in Table 1. By tracking the wall motion on both the anterior and posterior side of the blood vessel, the diameter distension curve can be determined as well as its associated distension measure ∆Dmax =Dmax-Dmin, with Dmax and Dmin respectively the maximal and minimal diameter during the cardiac cycle (cfr. fig. 5). Assuming planar deformation, radial (rr) and circumferential strain ( ) can be derived from the ultrasonic distension estimation ∆D as: D D rr     ) ( 